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Abstract 

We present a solution to the conservation form (Eulerian form) of the 
quantum hydrodynamic equations which arise in chemical dynamics by 
implementing a mixed/discontinuous Galerkin (MDG) finite element nu- 
merical scheme. We show that this methodology is stable, showing good 
accuracy and a remarkable scale invariance in its solution space. In addi- 
tion the MDG method is robust, adapting well to various initial-boundary 
value problems of particular significance in a range of physical and chemi- 
cal applications. We further show explicitly how to recover the Lagrangian 
frame (or pathline) solutions. 
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1 Introduction 



Quantum hydrodynamics (QHD) has engendered substantial activity in the field 
of theoretical chemical dynamics, where one may refer to Wyatt et al.([38 ) for 
a comprehensive introductory overview of the numerous recent results emerging 
from this blossoming field. 

The basic idea emerging from quantum chemistry in the context of QHD 
is to employ the time-dependent Schrodinger equation (TDSE) to solve for the 
dynamical properties (probability densities, "particle" velocities, etc.) of chem- 
ical systems. In the same spirit in which the de Broglie-Bohm interpretation 
(see [ll[5l[T6]) of quantum mechanics may be used to recover "trajectories" of 
individual fluid elements along the characteristics of motion of the solution, the 
QHD equations of Madelung and Bohm are derived as formally equivalent to 
the TDSE and thus comprise an alternative route to solutions which generate 
quantum trajectories that follow particles along their respective paths (see [3S] 
and [18j for a comprehensive overview). 




Figure 1: Here we have the intramolecular rearrangement of the aryl radical 
2, 4, 6-tri-fert-butylephenyl to 3, 5-di-fert-butylneophyl (see [6] for details). 

These solutions hold particular significance, where, in the context of the 
QHD formulation, it is possible to resolve the chemical dynamics of a vast 
number of reaction mechanisms known to have pathways dominated by quantum 



§1 Introduction 



3 




Figure 2: Here we show an enzymatic catalysis - an aromatic amine dehydro- 
genase (AADH) with a tryptophan tryptophyl quinone (TTQ) prosthetic group 
catalyzing the oxidative deamination of tryptamine with an electron transfer to 
an arsenate reductase enzyme (see [B] and [2^ for details, PDB codes: Inwp 
(azurin), 2agy (AADH)). 

tunneling regimes. Some of these systems include proton transfer reactions (for 
example see figure [ij, conformational inversions, biologically important redox 
reactions in enzymatic catalysis reactions (see figure [2]) , and proton-coupled 
electron transfer reactions (refer to j27j and |26]). It is not yet clear if these 
types of methods may also have application at higher energies, for example in 
the halo nuclei tunneling occurring in fusion reactions (as seen, for example, in 

Substantial research has been done in quantum hydrodynamics to find the 
best and fastest computational methodology for solving this system of equations. 
In the standard methodology presented using the quantum trajectory method 
(QTM), for example, solutions to the QHD equations are found by transforming 
the system of equations, which is generally posited in the Eulerian fixed coor- 
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dinate framework (see [T31 [TH [THl [H]), into the same set of equations in the 
Lagrangian coordinate framework, which effectively follows solutions along par- 
ticle trajectories; or along so-called "Bohmian trajectories." The transformation 
from the Eulerian to the Lagrangian frame leads to a set of coupled equations 
which solve for two unknowns: the quantum action S{t^ r) and the probability 
density or quantum amplitude \/ g{t, f) = R{t, r) along the trajectories r(i, a;) 
(e.g. see [SS] box 1.2). The obvious advantage of the Lagrangian framework 
is reduced computational times, since solutions are only computed along a set 
of chosen trajectories; while clearly the disadvantage is the possibility of ob- 
scuring structure hidden within the continuum of the full solution, which may 
only emerge properly in convergent numerical schemes, and also the increased 
complications of transposing into more complicated settings: such as with func- 
tional or time dependencies on the potential term V, or including dissipative or 
rotational vector fields. 

In addition, the numerical solutions to the above mentioned Lagrangian for- 
mulations have demonstrated characteristic behaviors which introduce certain 
technical difficulties at the level of formal analysis. First, the system of equa- 
tions are stijf, which is to say, solutions to the system may locally or globally 
vary rapidly enough to become numerical unstable without reducing numer- 
ically to extremely small timesteps. Furthermore, there exists the so-called 
"node problem," which is characterized by singularity formation (see [3S] for 
characterization of node types) along particle trajectories. Another issue which 
arises is obtaining unique solutions, since there is not a unique choice of tra- 
jectories in the Lagrangian formulation (see for example §6 and appendix A). 
And finally, boundary data is often treated without regard to the (often sub- 
stantial) numerical residuals introduced in the weak entropy case, or taking into 
account consistency between the TDSE and the QHD system of equations (see 
for example [H] and §3). 

We introduce an alternative formulation to the standard solutions described 
above in g and S and tracked with respect to the Lagrangian coordinate frame 
which is motivated by work of Gardner, Cockburn, et al. (see [TJ [TH [TS]). In- 
stead, we keep the system in its conservation form (instead of in a primitive 
variable form) in the Eulerian coordinate system (see |21]), and solve for the 
density g — g(t, x) and the particle velocity v — v(t, x) (instead of the quantum 
action S) . We show that these solutions may be used to easily recover the vari- 
ables S and ip in a, single step; and may with little difficulty be transformed into 
their Lagrangian coordinate frame counterpart solutions g{t,r),v(t,f), S(t,r) 
and 'ijj{t, f), using the conservation equation (continuity equation), or by solving 
for pathlines in the sense of classical mechanics, or by any number of alternative 
so-called "offset methods." Additionally, our solutions demonstrate a type of 
resolution invariance, which is to say that the behavior of our solutions are qual- 
itatively equivalent at varying spatial resolutions, and compare favorably with 
solutions to the formally equivalent TDSE. As a consequence, our conservation- 
based formulation is computationally competitive with Lagrangian formulations, 
up to a type of "formal accuracy" in the trajectory solutions. 

Our solutions, as the Lagrangian formulated solutions mentioned above, still 
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demonstrate a stiff behavior. However, also as the Lagrangian solutions above, 
and similarly to the classical CFL condition in fluid mechanics, we consider this 
a prohibitive but not insurmountable computational difliculty. On the other 
hand, our solutions to the conservation form of QHD do not demonstrate the 
node problem (at least on Gaussian wavepackets) as expected, as the only type 
of node our formulation exhibits is for g = 0, which never occurs if we add 
a numerical ambient density to the initial density p\t=o- The solution is 
stable when the ambient density is set to '--^ 11 orders of magnitude smaller 
than TaaxQ{g) over a computational domain fl. We maintain that the addition 
of to the initial density does not significantly change the numerical solution 
of the system of partial differential equations, while introducing the substantial 
benefit of significantly improving its stability. Again, this behavior compares 
favorably with solutions to the TDSE, which also do not demonstrate the node 
problem. On the other hand, computing solutions in the Lagrangian frame 
still offers substantial computational efficiency when compared to those in the 
Eulerian frame; due simply to relative density of solutions. 

We begin in §2 by presenting the governing equations, then rescaling these 
equations in time for substantial improvement of numerical tractability. Next we 
present the details of a computationally well-posed finite element discretization 
scheme leading to our approximate (numerical) solution. The scheme is based 
on a discontinuous Galerkin method for the QHD conservation laws and a mixed 
finite element method for the Bohmian quantum potential, which is inspired by 
[5]. In §3 we briefly derive the basic equations, and discuss the rather strong 
dependence on the formal and numerical equivalencies in the boundary data. In 
§4 we derive an analytic test case which allows us to find the relative error in 
the discontinuous Galerkin mixed method, which shows that our formulation is 
near to numerically exact everywhere but at the boundaries (which is expected). 
We proceed in §5 by testing the standard case of a hydrogen atom tunneling 
through an Eckart potential barrier, compare these results to a finite difference 
scheme for the TDSE, and then show how to use the continuity equation to 
recover the Lagrangian, or Bohmian, trajectories. Next, in §6, we show how to 
compute pathlines, recover the variables p, u, ip and S in both the Eulerian and 
Lagrangian frames, and compare the way in which these solutions relate to each 
other. We finish with some concluding remarks in §7. 

§2 Conservation Formulation of Quantum 
Hydrodynamics 

Consider the following system of equations for {s,x) G Tg x ft, motivated by 
|38j . where we have transformed the solution space from the usual Lagrangian 
coordinate frame into the conservation form of the Eulerian coordinate frame: 



dsigmv) + Vj, • n + gV^V = 0, 



(2.1) 
(2.2) 
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with initial conditions 

gs=o = go, and v^^o vq 

where g = g{s, x) is the probability density corresponding to conservation equa- 
tion (2.1 1, and v — v(s,x) is the volume velocity corresponding to the momen- 
tum density gp = gmv in equation (2.2 1, where the mass m is constant. Here 
V corresponds to the potential surface, where in keeping with the usual formu- 
lation in chemical applications in one dimension V may be generally thought of 
as a model potential (e.g. an Eckart, Lennard- Jones or electrostatic potential). 
The quantum stress 11 is given to obey, 



n 



gmv iSi V + g 



\4m 



4m 



Wig, 



or alternatively 



m = gv iS) V 



gh^ 
4m2 



log g, 



with the Bohmian quantum potential given as Q = ^^A^jy^^ / ^fg (note that 

this term is only defined up to a sign convention, see for example Ref. |18l I23j 
versus Ref. [38,), such that the nonlinear dispersion relation is given by. 



gt^ 
2m 



Am 



V-{gVllog g), 



yielding the alternative form of (2.2 1: 



dt{gmv) + • {gmv v) — gVxQ + g^xV = 0. 



(2.3) 



(2.4) 



Let us rescale (2.1 1 and (2.4 1 by setting s = y/mt and solving for a rescaled 
solution u and p in the time variable t, such that u(t,x) = ^/mv(y/mt, x) and 



(2.5) 
(2.6) 



p{t,x) = g{y/mt,x) such that (2.1) and (2.4) for {t,x) eT xfl become: 
dtp + "^x ■ (pu) = 0, 

dtipu) + Vx ■ {pu ®u)- pVxQ + pVxV = 0. 



We solve (2.5 1- (2. 6) using a mixed discontinuous Galerkin finite element 
method. We define the state vector 

U = {p,puf, 

the inviscid flux vector 



/ = {pu,pu (Ki u)^ 



and the source vector 



s^io,pVx{v^Q)y 

Then we can rewrite (2.1 )-( [2!2| as 

Ut + fx + S = 0. 



(2.7) 
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Figure 3: The discretization of f2, distinguishing nodes, elements and neighbors, 
with boundary dfl = {a, b} in dimension iV = 1. 



Consider the following discretization scheme motivated by [TH HH] (and il- 
lustrated in the one dimensional case in Figure |3|. Take an open f2 C M with 
boundary — F, given T > such that Qt = ((0,r) x fi) for fl the closure 
of il. Let denote the partition of the closure fi, such that taking f2 — [a, b] 
provides the partition 

a = Xq < Xi . . . < Xne = b 

comprised of elements Gi — {xi^i,Xi) G ^ such that = {Gi,G2, ■ ■ ■ ,Gne}- 
The mesh diameter h is given by h — sup^g^^ (a;^ — Xi-i) such that a discrete 
approximation to is given by the set fi^j = UiGi \ {a, b}. Each element of 
the partition has a boundary set given by dGi = {xi^i,Xi}, where elements 
sharing a boundary point dGi H dGj 7^ are characterized as neighbors and 
generate the set ICij = dGi H dGj of interfaces between neighboring elements. 
The boundary dfl = {a, 6} is characterized in the mesh as d^l = {xQ,Xne\ 
and indexed by elements Bj g dO. such that f2 = ^ U JCij U dO.. Now for 
/ C Z+ = {1, 2, . . .} define the indexing set r{i) = {j £ I : Gj is a neighbor of 
Gi}, and for Ib C Z~ = {—1, —2, . . .} define s{i) = {j ^ Ib ■ Gi contains Bj}. 
Then for Si — r{i) U s(«), we have dGi = ^j&s(i)^ij ^^nd dGi H dO, = yjj^s(i)^ij- 
We define the broken Sobolev space over the partition S/t as 

W'^-^inh, ^n) = {V : V\g^ e W'''\G^) yG^ E ■%} . 

Further, approximate solutions to (2.1|-(2.2| will exist in the space of discon- 
tinuous piecewise polynomial functions over restricted to S^u, given as 

Sti^h, ^h) = {V : V\g^ G 0^\G^) yG^ E ■%} 

for 3^'^{Gi) the space of degree < d polynomials on Gi- 

Choosing a set of degree d polynomial basis functions Ni G 3^'^{Gi) for 
£ = 0, . . . , d we can denote the state vector at the time t over fi/i, by 

d 

Uh{t,x)=Y,U\{t)Nl{x), Va;Ge„ (2.8) 

f=0 
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where the A^"s are the finite element shape functions in the DG setting, and the 
f/^'s correspond to the nodal unknowns. We characterize the finite dimensional 
test functions 

d 

ip„€W^'\nH,.%), by ¥'h(^) = E^^^^(^) 

where if\ are the nodal values of the test functions in each Qi. 

Assuming that the source term S is sufficiently smooth, we let [/ be a 
classical solution to (2.7) and multiply through by and integrating such 
that: 



d 
dt 



JGi JGi JGi 



(2.9) 



Integrating (2.9) by parts gives 
d 



dt 



f Uh-Vhdx+ f {f-Vf,)^dx-( f-ip''Jx^~~( S-ip^dx. (2.10) 



Let ipitCij and V'iKji denote the values of (p on ICij considered from the interior 
and the exterior of Gi, respectively. It should be noted that for ICij € T, the 
restricted functions cpf^\ICji are determined up to a choice of boundary condition, 
which we will discuss in more detail in §3. We approximate the first term in 
(2l0| by. 



d 
dt 



the second term using an inviscid numerical fiux by 



(2.11) 



N 

/ ^{fh)i- {n^j)l'Ph\K.,,dlC, 



(2.12) 



*j 1=1 



for Uij the unit outward pointing normal and where I is the dimension, and the 
third term on the left in (2.10) by: 



dx. 



(2.13) 



Using (2.11 l-(2.13l, taking the convention that 



and setting the inner product 



{ah,bh)ng= ^ / ah-bhdx, 
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we define an approximate solution to (2.9l-(2.17l as Uh for all t E (0,T) satis- 
fying: 

Discontinuous Galerkin Method for the QHD Conservation Laws 



1) UheC"{[0,T];S^^), 
3) Uh{0) = U„. 



(2.14) 



To compute the source term S, we approximate the Bohmian quantum po- 
tential using a mixed finite element method. In particular, we know that at 
each time t, the quantum potential Q satisfies the equations: 



Q 



2m yfp 



and q = V^^/jo. 



(2.15) 



Let d e L^(ri) and ^ G -ff(div,ri). Then multiplying (2.151 by i? and respec- 
tively, and integrating by parts over results in: 



Qddx ■ 



q ■ (,dx 



-ddx. 



n 2m 

^/pVx^dx 



fp^ ■ ndT. 



(2.16) 
(2.17) 



Choosing finite dimensional subspaces Jf'' C L'^{n) and J^'^ C i?(div, ft), a 
mixed finite element method for the Bohmian quantum potential is then: find 
Q'' : [0, T] X 17 ^ M, q'' : [0, T] x ^ ^^^-^ ^Yi&t for all t G [0, T], Q''{t) G 
and q'' G satisfy: 



Mixed Method for the Bohmian Quantum Potential 

1) iQ.,^.)n = '^(^,i 

2) {qh^'^h)n = -(VPft^, Vj;^h)o + (\/P^>?/i'^)i 



(2.18) 



Since we wish S G L^(ri), we choose to be a continuous finite element space, 
and we choose Jf'^ to be an i/(div)-conforming space (e.g. Raviart-Thomas 
elements |33j . such that in one dimension, Raviart-Thomas elements collapse to 



be standard continuous finite elements). Equations (2.141 and (2.181 define our 



mixed/discontinuous Galerkin method in semi-discrete form. Computationally, 
we must also discretize time, as shown in §4 and §5. 

It is worth noting that in the Lagrangian formulation the primitive variables 
(jO, u) are accompanied by the quantum action S and the quantum wave func- 
tion tp. We will explicitly derive these terms in section §5 from the solution 



(2.14). It is also worth noting that a pure discontinuous Galerkin method was 



implemented as an alternative approach to the MDG method solution shown 
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in (2.14). This treatment used a dispersive flux formulation as shown in 
We found that this formulation depended nonlinearly on the sign of the advec- 
tive flux term, leading in the naive implementation to the formation of soli- 
ton/compacton type behavior; solutions which are well-known in the 'formally' 
equivalent formulation of Korteweg fluids (see [31 [TTJ [Ml [H]) - up to turbu- 
lence effects etc., as explained in §3 - which model diffuse fluid interfaces as 
well as having phenomenological interpretation in the context of the nonlinear 
Schrodinger equation (see [35] ) and the Gross-Pitaevskii equation (see [H [19] ) 
given nearly identical initial conditions to the ones we use in §5. However, in 
the context of chemical dynamics it is not clear that these types of solutions 
carry physical significance, and so we have isolated our analysis to the MDG 
method formulation presented in (2.141. 



§3 Boundary Treatment 

A recurring difficulty in constructing numerical methods for initial-boundary 
value systems of partial differential equations for physical systems is the issue 
of how to prescribe mathematically consistent boundary conditions which ac- 
commodate dynamic (physical) boundary data. It turns out that this issue is 
a cause of both numerical and mathematical difficulties in establishing the for- 
mal equivalencies between the TDSE and the QHD system of equations. We 
show this behavior explicitly in an example in §5, but let us flrst examine the 
mathematical source of this difficulty. 



Recall that the system presented in (2.1 )-(2.2 1 is derived explicitly from the 
TDSE. That is, we have set ip = -Re*'^^ and want to expand the solution 
of the Schrodinger equation in one unknown and one equation in ijj = tpit^x) 
into a system of partial differential equations in the unknowns R — R(t, x) and 
S = S{t,x). To make this a well-posed system we of course need a system 
of two equations, where both unknowns must be assigned distinct boundary 
conditions. First take the following form of the Schrodinger equation: 

A,+ -^vU^—dttP, (3.1) 

and plug in -0 = Re^^/^ such that expanding gives for the time derivative. 



n 



^-II^e^s^%R~^-^Re'^^%S, 



(3.2) 



and for the spatial component 

= e'^/'' (a^R + |v,^V,i? - ^(V.^)2 + ^i?A,5 



§3 Boundary Treatment 



11 



Putting (3.2) and (3.3) back into (3.1) and canceling a factor of e^^^^ we 
obtain: 

R^-^V = A^R + ^-^rV^SV^R ~ §{^.Sr + '-RA^S + - 

(3.4) 



Now, collecting the imaginary parts of (3.4) 
2mi 



. dtR~'^V^SV^R~ '-RA^S ^0, 
n n h 

and multiplying through by h? /2m provides: 

-dtR - —S/^RS/^S - -^i?A^S' = 0. 
TO 2m 

Additionally multiplying through by ~2mR gives, 

m5ti?2 + V^R'^V^S + R^A^S = 0, 
where applying the product rule yields the conservation form: 

mdtR^ + V^ ■ (R^V^S) = 0. 



(3.5) 



Clearly setting R = y/g and using the Madelung relation v = ^'^xS for m a 
constant m e M leads to the usual conservation of mass equation: 



dtg + Vx ■ (gv) = 0. 



(3.6) 



Similarly putting together the real parts of (3.4) gives: 



2mR 



R 



^2 dtS~A,R+-{v,sy 



2m 
If 



RV = 0, 



such that upon multiplication through by h? /2m^R we have: 



m * 2m^R 



AxR ' 



2m 



TO 



Taking a derivation in x then yields 



1 



TO 



( 



\2m'^R 



AxR + V:, 



/ 1 



V 2to^ 



(V,5)- 



Now again we substitute the important Madelung relation v = giving 
the form: 

dtv + iv,(-u • v) - -i!-V,(i?-iA,i?) + -V,y - 0. (3.7) 
2 2m^ TO 

The Madelung relation, v = ^VajS", is of course equivalent to setting v to be an 
irrotational field, since for any field S, Vx x VxS = 0. Thus for an irrotational 
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vector field v, using that Vxiv ■ v) = 2{{v ■ Vx)v + v x Vx 'x v), we may rewrite 



(3.7 1 as, 

dtv + {v ■ Vx)v - —^Vx{R-'AxR) + -VxV = 0, 
so that multiplying by gm yields, 

gdtmv + {gmv ■ '^x)v - g—\/x{R''^AxR) + g\7xV = 0. 

2m 



Combining this equation with ( |3.6| yields: 

dtigmv) + \/x ■ [gmv ®v)- g^xQ+ Q^xV = 0, (3.8) 

for Q the Bohmian quantum potential given as Q = (]^I^x\/q) I \fQ- It is 
important to see that the formal equivalence between |3.1| and |3 8| is entirely 
dependent on Madelung's irrotational condition, which makes turbulent effects, 
for example, vanish. In the alternative derivation of the QHD regime, using 
moment expansions (see for example [T8', '38]) this restriction is not necessary. 
Thus we have arrived at our system of quantum hydrodynamic equations: 

dtg + ^x-{Qv)^^, ^gg^ 
dtigmv) + Vx ■ (gmv (^v) - gVxQ + g'^xV = 0, 

requiring initial conditions 

P\t=() = Po and U|t^o = "o, 

and numerically requiring explicit boundary conditions pf, and Ub on an irro- 
tational vector field v. Additionally, and as an important aside, the formal 
equivalence we have derived is constructed without mention of boundary con- 
ditions, which is satisfied over (0,T) x R^, but on a discrete domain C M'^ is 
a bit over optimistic, and as we will see below, does not in general hold. 

That is, the TDSE code (see §5) sets the initial data tpi^i, on the boundary as 
a time-independent condition, so the boundary value iph = ipb = tpi^b is enforced 
for all t G [0,r). Since t/^i ^ must be decomposed into Ri b and Si^b to make 



sense for the QHD formulation (3.9 1, these give Dirichlet conditions which can 
be implemented, but are unstable in the QHD regime, since Ri^b exponentially 
decays on the boundary and as a consequence is not numerically invertible; as 
it must be in the QHD formulation. These may however be approximated by 
setting pi^b — Pa, the ambient density, and Ui^b = ^mlg^^^'^^ the 
boundary element. 

However, these BCs still are not well-posed in the QHD regime for the fol- 
lowing reason. First we compute the entropy inequality for the rescaled ver- 



sion of n3.9\ shown in (2.5) and (2.6). We may compute the important classi- 



cal/quantum entropy satisfying for non-boundary terms that: 
d f f bp h^iVxJp? \ 
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We arrive at this system by multiplying the momentum equation from (3.9 1 
by V and integrating in space (e.g. the domain is some C M.^), such that 
rearranging we find 

/ vdt{pv) + vV X ■ {pv ® v)dx — I pvV xQdx + I pvVxVdx — 0. (3.11) 
Jn Jn Jn 

The product rule allows us to expand the first term on the LHS as: 
\v\'^{dtp+\7x ■ (pv)) + pvdtv + p\v\^^x ■ V dx, 



where |wp — v ■ v. Using the mass conservation equation twice from (3.9 1 and 
applying the divergence theorem we find that, 



V (dtipv) + V • {pv (g) v)) dx 



d_ 
di 



[ p\^dx+l [ \/x-ipv^)dx. (3.12) 
Jn ^ Jn 



Next, using the dispersion relation from the Bohm quantum potential the third 
term on the left yields: 



pv\/xQ = 



2m 
fi2 



pvVx 

1 



^xy/P 
VP 



dx 



/ Vj; • (y/pvAx^/p) dx 
Jn 



2m 



n 

2m 



Vx ■ {pv) Vx^/p+yx ■ i^/pvAxy/p) >dx 



2m 



■ [ \7x i^^x ■ {pv)^xVp) dx 
Jn \vP / 

/ ^ x\fp d-tSI x\fpdx + boundary terms. 
Jn 

Finally the source term Y — V(x) upon integrating by parts gives 

/ pv ■ VxVdx = - / Wx ■ {pv)dx + i Vx- {Vpu)dx 
Jn Jn Jn 

/ pVdx + Vx- {Vpv)dx. 
Jn Jn 



d 
dt 



(3.13) 



(3.14) 



Then we have recovered (3.101 as an equality up to the boundary terms in 



(3.12), (3.131 and (3.14). To recover the mathematical well-posedness of the 
system these boundary terms must either vanish or be bounded and positive (or 
negative) definite. One such choice of boundary data is, for example, Vi, = 0. 
Another is the pair of conditions Wx^/Pb = and Vf, = for all t E [0,T), and 
so forth. 

The first set of boundary data, with vi, = 0, may be set with pf, = pA- Since 
the action behaves as a phase, this seems a reasonable approximation, since it 
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effectively assumes that up to a constant of integration that the phase is constant 
over boundary elements V^b = 0. These conditions are then mathematically 



consistent with the system of equations (3.9), but have the physical effect of 
generating ^Hnlet/ outlet" boundary layers, caused by the value of pb. 
Perhaps a more natural boundary condition is given by setting, 

where C/JJ is the numerical solution at timestep as explained in detail in 
§4, and ICij G dfl. This boundary type is a first order approximation to a 
transmissive or radiative condition that treats the boundary like a "ghost cell," 
and allows density and momentum to leave the domain as though falling into 
vacuum, while allowing no density or momentum to enter. This condition ap- 
proximates to the first order, the effect of "not setting boundary conditions 



at all," and thus not badly perturbing the system (3.9 1 away from its natu- 
ral behavior, nor generating reflecting behavior, which in some contexts - such 
as a chemical reaction occurring in a solvent bath - are difficult to physically 
interpret. 

§4 A Numerical Test Case 

We wish to test the accuracy of our MDG method formulation by solving an 



analytic test solution. In order to do this we choose a numerical flux for (2.141 
and restrict to spatial dimension / = 1. For the inviscid flux $ we implement 
the local Lax-Friedrich's flux €>;lf satisfying 



/ ^iLF-'PhdlC^^ j {f{Uh)\K:,^+f{Uh)\K.,,)-nijcp,^\^^^dlC 

(Spec^(ro))((J7,0|K;.^ - {Uh)\jc,J ■ nij^h\K,jd^, 



1 

for riij the outward unit normal and Spec^(ro) the spectral radius of Fq; the 
Jacobian matrix of the inviscid flux Ju fiU) — Tq{U) which may be represented 
by the following 2x2 matrix. 

Summing over the elements of the mesh this term satisfies: 

die 



I (Spec,(ro))((C/0|;c.,-(t//,.)|^,J-n.,<P/.k.,rf/C. 



(4.2) 
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-7 Relative Error in fJ -s Relative Error on be 30 
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Figure 4: Here we show the relative error introduced by the weak entropy bound- 
ary conditions for a = and both 6 = 10 and b = 50. The boundary data (the 
graphs on the right) show only the relative error on element b of dCl for b — 10 
and 6 = 50, respectively. 



Next we discretize in time. That is, we denote a partition of [0,T] by 



for a timestep given as Ai" = — t' 
timestep t". Thus we implement the following forward Euler scheme: 



and let C/JJ denote the solution at 



Ul 



dt Ai" 

which, along with the implementation of a slope limiter in the conservation 
variables {p,pu) given by van Leer's MUSCL scheme (as shown in [36] and 
|57|). allows us to explicitly solve (2.141. That is, we define an approximate 
solution as for all t„ such that n — 0, . . . ,T satisfying: 



1) Ul e St Ql e and ql e 



2) 



h 



At" 



0, 



3) {Ql,^hh = 



5) U'i = Uh{0). 



(4.3) 
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The above formulation lends itself naturally to a staggered scheme. First, given 
f/JJ one solves step 3 and 4 for QJJ and q]^, which provides SJJ, allowing us to 
solve for U^~^^ in step 2. 

Now we construct an appropriate test case. Consider the dimension = 1 



case and let m = on for (2.5)-(2.6l, such that dsP = 0. Up to a choice of 



boundary conditions, upon integration we have for (2.6) that 

such that choosing a C = we find the following second order ordinary differ- 
ential equation: 

p"-p-i(pr-o, 



whose solution is p — e^. We solve for the approximate solution of (2.14) using 
the above scheme, with initial conditions po = e^,u = 0, V — C and m — 1836 
the mass of a proton in Hartree atomic units (au). The boundaries are set to 
the weak entropy boundary condition formulation as presented in [2 |3J US] and 
|29j . We graph the relative error of our approximate solution ph to the exact 
numerical representation pa in Figure ^. We see that the two solutions are 
numerically exact in the interior of the domain, and error accumulates in the 
boundary dfl, as expected due to the weak entropy boundary conditions. We 
note that the error on the boundary may be reduced by increasing the absolute 
size of the interval [a,b]. 



§5 Tunneling in TDSE and QHD 

We proceed by testing a relatively standard example in quantum chemistry, 
given by a propagating Gaussian packet in the direction of a model Eckart 
potential barrier. We solve the following one dimensional system: 

dtp + dM = o, (5.1) 

dtipu) + d^ipu^) -pdM + pd^V = 0. (5.2) 
with initial conditions 

Po^ PA+ { Ir — ) e 2^"°' and uq = (aVb)^^^ , (5.3) 



/27r/i 

where the Eckart potential is given by 

V{x) = Vo sech^ Q(a; - xi)j . (5.4) 

As is conventional in quantum hydrodynamics, the mass is set to approx- 
imate the hydrogen (proton) mass m ~ 2000 au (in Hartree atomic units), 
PA ^ 10^^" is a numerical background density for division, xq centers the Gaus- 
sian packet, xi centers the potential, p is the variance of the distribution, a 
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TDSE Finite Difference Solution QHD Finite Element Solution 




Difference in the TDSE and QHD Solution 




400 



Figure 5: The top graphs compare solutions to the TDSE and QHD system in 
the so-called "eyeball norm," for the forward Euler scheme. The bottom solution 
shows the nontrivial formal difference. Here x refers to the a;-th meshpoint. 

is a constant a G M and Vq is the barrier height (which we may vary, so some 
constant Vq € K). In the quantum regime (when classical barrier transmission 
is not present), the initial velocity uq is often chosen to satisfy the following 
condition on the initial kinetic energy Kq = }^u^= jVq. 

The background ambient value pA is required in order to satisfied the math- 
ematical and numerical well-posedness of the system such that the behavior of 
the system is not perturbed away from its proper character by compounding 
residual behavior, as shown in §3. Furthermore, from a phenomenological point 
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min(p,10-=) in QHD min(p^„p,10-=) - min(p_^,10-=) 




400 800 1200 1600 2000 400 800 1200 1600 2000 



f(au) ;(au) 

Figure 6: We show the diffusive noise profile niin(pQHD, 10^"^) in the QHD 
solution, and the difference min(pxDSEi 10^^) — min(pQHDj 10^'^). Here x refers 
to the x-th nieshpoint. 



of view, this value is nonrestrictive and physically easily justified - for example, 
for a chemical reaction occurring in a solvent bath, or, similarly, any process 
occurring away from vacuum. 

The discretization proceeds as in section §2 and §4, where we adopt the local 
Lax-Friedrich's inviscid flux with van Leer's MUSCL slope limiting scheme. 
Next we implement a standard explicit Runge-Kutta time discretization (see 
[21 [34] and [29], or [28] for explicit details). 



Now we solve the resultant system using for our initial data (5.3l-(5.4) ex- 
plicitly that fi = 0.16, a — 2,xa — 3 and xi = 6, such that. 



po = 10-1" + ( — ^ ) and uq = (2Vb)^/^ 



with potential: 



V{x) = Vq sech^ 



It is worth noting that we have thus chosen a kinetic energy which is in the 
context of a mixed classical-quantum regime; which is just to say that some 
classical trajectories trasmit over the barrier, in addition to those that tunnel 
quantum mechanically. For boundary data we use the approximate well-posed 
Dirichlet conditions discussed in §3: 

Pb = PA = lO^i^" and Uf, — 0. 

We compare our solution to a finite difference scheme for the TDSE provided 
by Prof. Robert E. Wyatt [38^ in order to test the accuracy of our formulation. 
The TDSE has equivalent initial settings, while the boundary conditions are 
given naturally via ^pt, — tpi^b as discussed in §3. 

In Figure [5] these two solutions are compared. Is is clear that the two solu- 
tions have the same qualitative behavior. However they do show fundamentally 
different quantitative behaviors. Analysis has shown that the two most preva- 
lent sources of error between these two solutions are diffusion and boundary 
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Absolute Difference in QHD Boundary formulations 
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Figure 7: We show the absolute difference between the QHD solution using 
the approximate boundary data from Figure [5] denoted with the transmis- 
sive boundary formulation from (5.8 1 denoted pr- Here x refers to the a;-th 
meshpoint. 



oscillations. The boundary oscillations clearly occur due to the approximations 
discussed in §3. The diffusion, on the other hand, is a signature of the slope 
limiter in the QHD formulation and is shown in greater detail in Figure |6] Here 
we confirm that the MUSCL slope limiting scheme is adding a type of "artificial 
diffusion" to the QHD solutions. We have found that choosing a less restrictive 
slope limiter, such as the flux limiter of Osher presented in [21], does stably 
reduce the diffusion in our solutions. 

We may now recover trajectories, or characteristics, of the solution by using 



the fact that (2.1 1 is satisfied at every time step (note that we show the alter- 
native method of integrating velocity "pathlines" in §6). We may think of this 
equation as a kind of "conservation of density" here, and thus we simply employ 
Reynold's transport theorem (RTT): 



d_ 
di 



pdx + / purei ■ ndx — 0, (5.5) 

f2(t) Jf(t) 



where Urei is the relative velocity of the fluid with respect to the moving bound- 
ary T{t). First consider the case when u{a) « such that we may choose 
fl(t) = {a,y(t)) where y(t) is the moving boundary treated as an unknown. By 
assumption and construction, Urei{a) = 0, whereas for a trajectory we require 



Urei{y) = 0. Then integrating (5.5 1 in t we flnd 



vW ry{o) 

pdx = / pdx. (5.6) 
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Accumulated Density Trajectories 
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Figure 8: We solve the accumulated density trajectories from (5.6 1 using the 
transmissive solutions from px in Figure |7] 



Let us define for each trajectory y{t) with y(0) = j/q the "locally accumulated 
mass" M by: 

fV{t) 

M{yo,t)^ / p{x,t)dx. 



Approximating each trajectory then directly follows from the equation Af (j/o, t) = 

M(yo,o). 

To continue let us denote Mi(t) — M(xi,t), where Xi is the i-th meshpoint. 
To compute y{t), we compare -M(yo, t) to the increasing sequence {-Mi(i)}i=o...Af 
and find j such that Mj_i(i) < M{yo,t) < Mj{t), which gives us that y{t) G 
). Then to find y{t) recall that we have from (|2.8| an expansion 



Ph{t,x) ^^ci{t)N^{x), for xe{xj-i,Xj) 

1=0 

where the c; = ci(t) are constants for every fixed t and the shape functions 
(x) in our implementation are translated versions of polynomials {Pi}f^Q on 
[—1,1]. That is using fj : [ ] I—!- [—1,1] where 

/,(x) = 2f^^^)-l, 

\Xj Xj^l / 

we find, (x) = Pi{fj{x)). Then solving for y{t), formulated via 

rvit) rv{t) d. 

M{yo,t) = Mj_i(t) + / Phix,t)dx = Mj_i{t) + / }2ciPi{fj{x))dx, 

can be recast by a change of variables, as solving for X in 



M(yo,t) = Af,_i(<) + 



2 



- \ I y,ciPi{z)dz, (5.7) 
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Figure 9: Here we show mass conservation in the QHD regime given transmissive 
boundaries, where bf = j,.. Jg^ pdxdt is the boundary flux. 



after substituting z = fj{x). But that just corresponds by the change of vari- 
ables, to 

Then a solution to X exists by the intermediate value theorem, and since the 
integrand is positive it is uniquely determined as the only solution on [—1, 1] to 



the polynomial equation of degree d + 1 arising from (5.7). We my then, for 
example, in the piecewise linear case (i.e. d = I) use the quadratic formula to 
recover X and hence the position of y{t) within Qj. 

Similarly, we also work in the other direction, with the balance of the mass 
in [y(t),Xj] such that the analogous integral equation becomes: 



/I d 



which provides for a consistency check on the accumulated density in either 
direction. Consequently we have that the sequence {y(i)}t=i,...,T provides a 
numerical approximation to the position of a particle initially at yo when t = 
at our given set of later times. 

This formulation holds as long as our hypothesis, u{a) « 0, is satisfied. 
However, we can immediately extend this result to include the case u{a) ^ 0. 



That is, after integrating in t we note that (5.5 1 becomes 



y{t) ry{0) rt 

pdx = / pdx + / pu{a)dt. 
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Figure 10: Here we show the remarkable spatial invariance of the solution. These 
represent the same solution as that given in Figure |5] except the left graph is 
with 25 meshpoints and 100 timesteps, and the right at 50 meshpoints and 200 
timesteps. 



This gives us an alternative equation to find y{t), where we must only add 
pu{a)dt to the accumulated density M at every t. We further note that this 
basic framework may also be adapted to higher dimensions (see |32j). 



Now, we again solve our system with (5.6 1 using for (5.3l-(5.4) that fi = 
0.16,0; — 2,xo = SjXi = 6, however now we introduce the transmissive bound- 
ary condition: 

U^;;\^^^ = Ul\^^^, (5.8) 

as discussed in §3. In Figure[7]it is clear that the behavior between the solutions 
with transmissive and approximate solutions is quite distinct, and that bound- 
aries are, so to speak, felt in the interior solution even before significant density 
has reached dfl. We use the transmissive boundary conditions to construct the 
accumulated mass trajectories derived above, as they seem to represent more 
physically cogent boundaries. The results are shown in Figure |8] where the 
the "Gaussian centered trajectories" are simply the trajectories containing the 
majority of the initial density; that is, those trajectories whose initial positions 
are ±1 au from the center of po- 

We further show that the MDG method is a conservative scheme. That is, 
in Figure [9] we show that the density is effectively (numerically) normalized to 
one on the domain, when taking into account the boundary flux. That is, we see 
linear error growth at machine precision over 10,000 timesteps. Another feature 
of the solution which is attractive in the sense of practical applications, is that 



the spatial invariance demonstrated by the solutions. In Figure 10 we show 
the this feature, where the same calculation from Figure |5] is graphed, where 
there 400 meshpoints and 10,000 timesteps were used in order to compare with 
the TDSE. However, as is clear from Figure [TO] with only 25 meshpoints the 
solution provides the same qualitative answer. This is an important feature in 
chemical applications where computations must scale in SN dimensions, for N 
the number of atoms in the molecular system of interest (see for example pO|). 
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§6 Recovering and S in both frames 

Now that we have solutions in the Eulerian and Lagrangian coordinate frames as 
given in §4 we may recover the important variables ip and S in either frame. First 
we note that we may alternatively recover the trajectories using the solution U 
from §4 to solve the initial value problem: 



dr 
dt 



= u{t,r) with r^^Q 



(6.1) 



We recover these r by direct integration, and compare them to those computed 



via (5.6 1 (see figure 111, where we refer to the r trajectories computed in (6.11 
as the "velocity pathlines." 
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Figure 11: We graph the quantum trajectories using (6.11 to solve for r, which 



can be compared with the accumulated mass trajectories shown in Figure |8] 



The trajectories computed using the velocity field (6.1 1 are shown in Figure 



|11| and show qualitatively similar behavior to the trajectories computed using 
the accumulated mass formulation in (5.6 1. There is no necessarily unique 



way of arriving at the trajectories one chooses to represent the solution in the 
Lagrangian frame. For example, one may utilize a method which weights the 
solutions between (5.6 1 and (6.11. That is, we may compute the trajectory 



positions via (6.11 and then offset these by a weighted average of the density 




Figure 12: The Eulerian solution p{t, x) and the corresponding Lagrangian so- 
lution p(t, r) for the same initial condition settings as in figure [S] using the 
conservation form of the trajectories (5.6 1.. 



24 



Chemical PDEs 



Eulerian Frame Soluton 

200- 
100 


-100 

1000 2000 

;(au) ^(au) 



Lagrangian Frame Soluton 




Figure 13: A graph of the Eulerian solution S{t,x) and the corresponding La- 
grangian solution S(t, f) for the same initial condition settings as in figure 12 



using the conservation form of the trajectories (5.6) 



conservation in (5.6 1. We provide details on particular alternative in appendix 



A and show an example case. 

It is now possible to solve for a number of derived variables in either the 
Lagrangian or Eulerian frames in order to recover the phase information of 
the quantum wave-packet associated to each characteristic pathline. First we 
recover the trajectory-wise solutions p{s,r) and u(s,r), and then compute the 
variables: 



VxS — \frnu and ■0 



I pe 



iS/h 



(6.2) 



where S{s,r) is the quantum action and ip{s,f) is the quantum wavefunction. 
Recall that R{s,'r) — \/p{s, r) as shown in §3, such that using v for the velocity 
from §2 we recover from (6.2) the more familiar formulation: 



and V = Re'^/^. 



(6.3) 



It is important to note that up to a constant of integration, S and ip are 



completely determined by the solution (2.14). Also, (6.3) is satisfied in both 



reference frames, so we now have the following solutions: 

p{s,x),p{s,r),v{s,x),v{s,r),tlj{s,x),ip{s,r),S{s,x), and 5(s,r). (6.4) 



These solutions are graphed in figures [12] - fr3| where it is interesting to note that 
the two frames draw out different aspects of the solution. While the Lagrangian 
frame tracks individual "particle" trajectories across the function profiles, it 
misses some of the nuance in the continuous structure of the surface; which 
is naturally recovered by the Eulerian frame solution. Furthermore, as lower 



resolution, we find that the conservation based trajectories from (5.6) are more 



well-behaved than the velocity based trajectories from (6.1) 



§7 Conclusion 

We have presented a numerical solution to the quantum hydrodynamic equations 
of motion as posited in the context of quantum hydrodynamics with chemical 
applications. Our approximate solution is a rescaled (in time) version of the 
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standard QHD equations and is the first model of its type presented in a mixed 
discontinous Galerkin framework in the context in which it arises in chemical 
applications. Our solution further shows good stability, up to a stiffness of 
the system of equations which is a well-known feature of the QHD system of 
equations, and a scale invariance behavior which makes it very appealing for 
the so-called "fast and dirty computations" often needed in realistic chemistry 
applications. Additionally we have shown in a rigorous and consistent way 
how to prescribe proper boundary data, which is often bypassed in the usual 
Lagrangian formulations of the system. We have further demonstrated that in 
the conservation formulation of this system, the quantum wavefunction ^ and 
quantum action S, which are used as motivation for the derivation of QHD 
systems to begin with (e.g. |51[211[3S]), are in fact completely determined (up 
to a constant of integration) by the solutions g and v. 

Finally, it is worth mentioning that these solutions are very closely related to 
quantum hydrodynamic solutions which have been extensively studied in other 
fields (see [ini [13 HH UHl 123]), but still maintain some important differences. 
One of the most important and prohibitive aspects of the quantum chemical 
formulation of QHD, is that the potential surface V arises from a multiple of 
3N degrees of freedom of each quantum subsystem, for N the number of atoms 
in each molecular subsystem (for example in an intramolecular rearrangement). 
This arises from the interpretation of the wavefunction ip as being the foun- 
dational variable in the dynamics of the quantum subsystem in the chemical 
models. Clearly, even for relatively small molecules, this immediately leads to 
extremely difficult numerical problems. In this sense it is important to have 
a numerical scheme which is easily parallelizable, fast, robust and accurately 
reflects the mathematical character of the solution. The MDG formulation pre- 
sented herein is a numerical method that fulfills these requirements, and offers 
a viable solution to some of the many difficulties which arise in the complicated 
solution space of chemical quantum hydrodynamics. The scale invariance of the 
solution makes it even an alternative approach to the Lagrangian formulation; 
up to the "formal" accuracy of solutions. 
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Appendix A 



The conservation method of recovering trajectories in (5.6 1 and the velocity 



integration method of recovering trajectories in (6.1 1 in no way exhaust the 
number of ways of representing solutions in the Lagrangian frame. In fact, 
there are an infinite number of ways of choosing trajectories. We introduce a 
way of computing a subset of these, and refer to these as "offset methods." 




Figure 14: On the top we show the quantum trajectories using the offset method 
solution of the same problem in Figure 11 with r = 1; and on the bottom we 
show the same trajectories using r — 2. 



The offset solution relies mainly on velocity integration but includes some 
information from mass conservation as follows: velocity integration provides an 
estimated position for each particle at the following time-step. Then one works 
through particle by particle, starting at the new estimated position and using 
mass conservation to estimate the new positions of its neighbors (a tunable 
number of consecutive elements on either side) offset from the velocity estimate 
of the 'current' particle. We set our tuning parameter to r here on both sides, 
though there is no reason a priori to choose a symmetric (with respect to ei- 
ther side) tuning parameter. Generically this provides a set of estimates for 
the position of each particle: one directly from integration, and others via the 
relationship of that estimated position to the relative estimated position of its 
nearest neighbors. 

That is, if P™ is the velocity estimated position, and P^-r P^+r ^^'^ 
the positions of the particles on either side that density conservation requires. 
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and applying our symmetry constraint gives for the new position that: 

r 



where the wiS are the weights for each component, in our examples computed 
with a Gaussian weighting functions uji — e~('"(^)/'' such that: 



Wi = LOi I LOi for i = 0, 

1=0 



Then for r = 1 we have = 1/2 and wi = 1/4. We show two examples of 
obtained offset trajectories in Figure [14] which are located at distinct locations 
in the solution space. Also note that these trajectories behave substantially 
different than those in Figures [8] and 1 1 
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